rm(list = ls())
options(warn=-1)
library(rgl)
library(misc3d)
knitr::knit_hooks$set(webgl = hook_webgl)

Load plot3D functions.

source("./plot3D_func.R")

STitch3D’s spatial domain detection result reveals layer structures in cortex.

celltypes <- c("Ext_Hpc_CA1", "Ext_Hpc_CA2", "Ext_Hpc_CA3", "Ext_Hpc_DG1")
celltype_colors <- c("#d62728", "#FFFF00", "#1f77b4", "#2ca02c")
um <- c(0.98223096, -0.05770408, -0.1785820, 0,
        0.02726252, -0.89759272, 0.4399813, 0,
        -0.18568322, -0.43703181, -0.8800704, 0,
        0.0000000, 0.0000000, 0.0000000, 1) #set the initial view of the 3D plot

plot3D_proportions(directory = "./res_mouse_brain",
                   celltypes = celltypes,
                   celltype_colors = celltype_colors,
                   um = um)

#3D coordinates and model are based on CCFv3.
load("spotstable.RData")
load("VOLUMESMALL.RData")
load("VOLUME.RData")
for (i in (0:34)){
  prop <- read.table(paste0("./res_mouse_brain/prop_slice", i, ".csv"), sep=",", header=TRUE)
  if (i == 0){
    prop_all <- prop
  }else{
    prop_all <- rbind(prop_all, prop)
  }
}
library(stringr)
prop_all$X <- str_split_fixed(prop_all$X, "-", 2)[, 1]
spots.table <- spots.table[prop_all$X, ]
spots.table$X <- rownames(spots.table)
spots.table <- merge(spots.table, prop_all, by=c("X"))

prop <- spots.table[, celltypes]
prop$max_prop <- apply(prop, 1, max)
prop$max_celltype <- apply(prop, 1, function(x){celltypes[which.max(x)]})

alpha_threshold = 0.2
um <- c(-0.6389505, -0.1887468, -0.7457324, 0,
        0.4022625, -0.9083022, -0.1147684, 0,
        -0.6556883, -0.3733116, 0.6562859, 0,
        0, 0, 0, 1)

open3d(windowRect = c(0, 0, 720, 720))

## glX 
##   3

par3d(persp)

## NULL

view3d(userMatrix = matrix(um, byrow=TRUE, nrow=4))
drawScene.rgl(list(VOLUMESMALL))

for (c in 1:length(celltypes)){
  idx_celltype = (as.vector(prop$max_celltype) == celltypes[c]) &
                 (prop$max_prop > alpha_threshold)
  spheres3d(spots.table[idx_celltype, ]$AP.paxTOallen - 530/2, 
           -spots.table[idx_celltype, ]$DV * 1000/25 - 320/2, 
           spots.table[idx_celltype, ]$ML * 1000/25, 
           col = celltype_colors[c], radius=5, 
           alpha=1)
}